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We study the head-on collision of black holes starting from unsymmetrized, Brill-Lindquist type 
data for black holes with non-vanishing initial linear momentum. Evolution of the initial data is 
carried out with the "close limit approximation," in which small initial separation and momentum 
are assumed, and second-order perturbation theory is used. We find agreement that is remarkably 
good, and that in some ways improves with increasing momentum. This work extends a previous 
study in which second order perturbation calculations were used for momentarily stationary initial 
data, and another study in which linearized perturbation theory was used for initially moving holes. 
In addition to supplying answers about the collisions, the present work has revealed several subtle 
points about the use of higher order perturbation theory, points that did not arise in the previous 
studies. These points include issues of normalization, and of comparison with numerical simulations, 
and will be important to subsequent applications of approximation methods for collisions. 

PL, " 

\o 

I. INTRODUCTION 

The accurate prediction of gravitational waveforms produced in the collisions of black holes has become a central 
topic of research in general relativity, due to their potential observability with modern interferometric gravitational 
wave detectors. Given the lack of symmetries in a collision, it was believed for a long time that only a full numerical 
integration of the Einstein equations would lead to reliable answers. Recently it has been noticed Jl],|| that one can 
qq ' make some progress in understanding the collisions by using black hole perturbation theory, especially for collisions in 
which the holes start sufficiently close to each other for the collision to be considered to be the evolution of a single, 
distorted black hole. Following this approach, called the "close limit approximation," linearized perturbation theory 
q has been shown to provide a remarkably accurate picture of the head-on collision of momentarily stationary 0J^] and 
i ■ boosted black holes || . 

^Jy If this technique is going to be considered a valid method for cases in which full numerical simulations are still 
not available, one needs to develop indicators for deciding when the approximation is trustworthy. Intuitive rules 
of thumb, such as requiring that a single, almost spherical, horizon initially surround both holes turn out to be 
, too conservative to be practical, as was demonstrated by the momentarily-stationary head-on collision results [Q. 
5^ 1 Recently 0, it was suggested that the use of second order perturbation theory to provide "error bars" could be an 
effective way of estimating the domain of validity of the first order results. For the head-on collision of momentarily 
stationary black holes the proposal appeared to work very well Q . 

The purpose of this paper is to explore the application of second order calculations to an important case of initial data 
that is not momentarily stationary, the head on collision of initially "boosted" (i.e., moving) holes. The introduction 
of boost turns out to add several technical and conceptual complications. These are important beyond their relevance 
to the specific collision studied here, since the cases of realistic physical interest (collisions with spin and net angular 
momentum) all involve initial data that is not momentarily stationary. T his p aper will therefore attempt to lay part 
of the groundwork for future investigations of more realistic situations [T^Jl9[ |. 

One of the conclusions of paper will be that perturbative calculations can be a very reliable tool to get quantitative 
predictions at a certain level of accuracy. If we wish to push the accuracy to a few percent level, some questions 
remain. This is in part due to the fact that numerical codes we use for comparison cannot at present be trusted to 
that level of accuracy either, and in part due to several technical complications that appear in perturbation theory. 
It is remarkable, however, that the addition of considerable amounts of boost to the black holes does not preclude the 
applicability of perturbation theory techniques. 

This paper is closely related to two previous studies to which reference will frequently be made. The first is the 
linearized analysis for initially boosted holes, by Baker et al. we shall refer to this as BAABRPS. The present 
work relies heavily on the second order formalism described in Refs. ( |]l4],[I|Jl3[ ) , which we shall refer to collectively 
as GNPP. 
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In the present paper, Sec. II gives the details of how our initial data is parameterized with a separation L between the 
holes, and a momentum parameter P. A discussion is given of the meaning of perturbation theory in a two parameter 
space of solutions. It is also pointed out that a feature of our initial data differs from that used in computations 
with numerical relativity, and this difference hinders a perfect comparison of results. We find, in this section, that 
to use the second order formalism of GNPP, we must make gauge transformations that eliminate the first order 
monopole perturbation of the extrinsic curvature. This is the first of several technical issues that were not evident in 
the linearized work of BAABRPS or the momentarily stationary data of GNPP. 

Section III shows how the initial data is evolved with a wave equation that has the structure of a Zerilli || equation 
with a source term quadratic in the first order perturbations. From these results, it is shown how one computes the 
gravitational waveforms and energies correct to second order, and the second-order correct results are presented and 
are compared with results from numerical relativity. In this section a discussion is given of a second order technical 
detail that has previously been ignored and that was of little consequence in the evolution of momentarily stationary 
initial data. The gauge fixing used in GNPP leaves unfixed a degree of freedom associated with time translations. This 
is not relevant to the computation of radiated energy, but must be resolved if our waveforms are to be compared with 
those of numerical relativity. In Sec. Ill a convenient method is given for fixing this gauge freedom in the waveforms. 
With this choice, the second order correct waveforms as well as energies are found to be in excellent agreement with 
the results of numerical relativity. 

The methods and results of the paper are briefly reviewed in Sec. IV and the connection to future work is pointed 
out. 

We use for the most part notation introduced in BAABRPS and GNPP, which in turn is based on the notation 
of Regge and Wheeler Eg ]. In addition, we will use here the convention of adding superscripts in parentheses when 
necessary to indicate whether a quantity is first or second order, and to indicate what multipole it refers to. Multipole 
indices will be distinguished with £ =. Thus, for example, would indicate the third order quadrupolc 

Regge- Wheeler perturbation H . 



II. INITIAL DATA 



A. The conformal approach 

The momentarily-stationary Misner JtJ initial solution to the initial value problem of general relativity for a two 
black hole situation has a convenient explicit analytical form; for initially moving holes no such form is available and 
the first step in the problem is to find an appropriate initial value solution. We use the conformal approach (see S and 
references therein), in which one assumes the metric to be conformally flat g a b = ^ 4 <5 a fc and constructs the conformal 
extrinsic curvature K ab = <j) 2 K ab . In terms of these variables the initial value constraint equations (assuming maximal 
slicing TrK — 0) read, 

v a k ab = o (i) 

1 K ab K i. 

(2) 

where all the derivatives are with respect to flat space. 

One can construct B solutions to the first set of equations (momentum constraint) for a single black hole with 
linear momentum, 

KT = ^[2P(an b) ~(S ab ~n a n b )P c n c ] . (3) 

Here r is the distance, in the conformally related flat space, from the origin and P a is a vector in that space and 
can be shown to be the momentum of the hole in the asymptotically flat physical space. By superposing two such 
solutions one obtains the conformally related extrinsic curvature representing two moving black holes (although the 
flat space vectors P a in this case must be considered parameters that have a clear interpretation as momentum only in 
the case that the holes are widely separated). Since the constraint equation ([[]) for the conformally related extrinsic 
curvature is linear, the superposition still solves the constraint. As was done in BAABRPS, we locate the two black 
holes on the z axis of the conformally related flat space, at positions z — ±L and we choose the flat space vectors P a 
to be symmetrically directed towards the origin and to have equal size P. (The case of holes moving away from the 
origin is represented with a negative value of P.) 



2 



As in BAABRPS, we treat the separation parameter L as our perturbation parameter, and we expand the two- hole 
superposed K a b in L. Due to the equal mass/opposite momentum symmetry, only odd powers of L appear, and the 
first two terms are: 
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Here R is the flat space distance to the origin, related to the flat space distance d to the holes by d = 
\J (R 2 ± 2R cos 9 + L 2 ) , and the expressions in Eq. (Q) are valid only for R > L. 



B. Boundary conditions 

One must now put this solution in the right hand side of the Hamiltonian constraint (||) and solve the resulting 
nonlinear elliptic equation. In this process one needs to decide which boundary conditions to impose for the elliptic 
problem. A common choice in numerical studies has been the use of symmetrization of the data through the two 
throats of the black holes (see for instance || and references therein). This kind of procedure is not very convenient 
if one is interested in semi-analytic work as we are, chiefly because symmetrizing implies using the method of images 
an infinite number of times and the expressions involved become quite large and difficult to handle. On the other 
hand, nothing prevents one from constructing unsymmetrized data for boosted black holes along the same lines as for 
the momentarily-stationary case (the Brill-Lindquist problem |Io|] ) . This was recently emphasized by Brugmann and 
Brandt pl| . Here we will take this latter approach. This is not completely inconsequential, since the only numerical 
simulations available for comparison are for symmetrized data; we will return to this point later. To generalize the 
Brill-Lindquist construction to the case with momentum, one assumes the conformal factor to be composed of two 
pieces, 

(f> = Acg + 4>BL- (5) 

One piece 



m m 



\R-Ri\ \R-Ra 



(6) 



is singular at the points R — R\,2 in flat space, and represents throats. That is, when one introduces a new radial 
coordinate of the form 1/\R — Ria\ the "singular point " R = R12 of the conformally related flat space is seen to 
have the actual geometry of space that is asymptotically flat as \R — Ri a| — ► 0. The result of putting Eqs. (^|) and 
is 

V reg = --j- —7 rf , 7 

8 (4>BL + <P Icg )' 

which is to be solved for <^ reg with the boundary conditions that 4> leg is regular at R = Ri^ and approaches unity as 
r — * 00. Notice that the right hand side of Eq. (0) is well behaved at R = R%,2] although the numerator diverges as 
\R — -R1.2I 3 , the denominator increases as \R — Ri,2\ 7 ■ The main difference between our approach and that of Brandt 
and Brugmann |ll|| is that we shall solve the initial value problem perturbatively. 



C. Perturbative issues 



We have defined a two parameter (P and L) family of initial data that we could evolve into a two parameter family 
of spacctimes. We are, of course, primarily interested in the close limit, the limit of small initial separation, and hence 
of small L. In principle, we could use initial data correct to second order in L. This would mean solving Eq. (^) for 
</> reg and expanding the solution in L. In practice, Eq. (M) would require numerical solution since the right hand side 
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of Eq. (|t]) is regular, and the Green function for the equation is simple, this would present no significant obstacle, but 
it would have the disadvantage that we would have the solution only numerically. In particular, this would mean that 
the dependence on P would not be transparent. For that reason, we follow a different path. As in BAABRPS, we 
consider small P and small L. More specifically, we consider a P — r\L curve in the family of spacetimes, with rj a 
numerical factor of order one. This means, for example, that terms proportional to P 2 , PL and L 2 are all of the same 
order, and are our lowest order perturbations. Our second order perturbations will be of the form P 4 , P 3 L, ■ ■ ■. Due 
to the symmetry of our configuration, no terms arise of order P 3 , P 2 L, ■ ■ ■. 

Since K a b has a leading factor of P, the numerator on the right hand side of (||) is proportional to P 2 . This point 
is rather subtle, and there is a temptation to come to a wrong conclusion. The expressions in (|4|) are O(PL) and 
suggest that the numerator in (||) is 0(P 2 L 2 ). It must be understood, however, that the expressions in (|J) are valid 
only for R > L. The solution of does not depend locally only on the large R form of the right hand side. It 
connects boundary conditions at infinity with boundary conditions of the throats, boundary conditions for which the 
R > L condition does not hold. Due to this non- locality in (J2J) , or equivalently (Q), the conformal factor depends in a 
complicated, nonpolynomial, way on the P parameter. Our "small P" assumption amounts to taking the numerator 
in (g), or equivalently (0), to be perturbative, at all points in space. This, and closely related issues, are further 
discussed in BAABRPS. 

The perturbative problem requires at several points the specification of a "mass" , either of the spacetime or of 
the black holes. Let us discuss this in some detail. First, there is the problem of what mass does one use for the 
background Schwarzschild spacetime around which we are doing perturbation theory. Our experience shows that one 
should use the ADM mass of the spacetime. We saw a similar situation when we analyzed the radiation generated 
as an initially conformally flat ( "Bowen-York" ||) spinning hole |l^] settles into its Kerr final state. The spin rate 
was taken to be small, and the problem was treated as a perturbation away from the Schwarzschild geometry. For 
apparently moderate amounts of spin, the radiation generated was rather small, but the effect on the ADM mass 
(i.e., the spin dependent increase over the Schwarzschild mass) could be a factor of several. By computing exactly 
the effect of spin on the ADM mass, we found we could successfully apply perturbation theory for moderate spin. The 
only question could be if one uses that of the initial slice or that after the radiation has gone out, but for the cases 
of interest the difference is less than 1%, so we will consider the ADM mass of the initial slice for our background. 

Then there is the issue of the initial data. The initial data for boosted black holes is characterized by the separation 
L, the momentum P and a "bare" mass for each hole to, which also serves as overall scale factor. This mass has no 
clear physical meaning, and no equivalent in the reflection symmetric initial data used in numerical relativity. Because 
of this, we would prefer to have the initial data parameterized by P, L, Madm- Since P, L and to determine uniquely 
the ADM mass, this is formally no problem. In practice we proceed in the following way. The ADM mass (for a given 
set of parameters L, P, to) can be found from the monopolc part of our second order solution for the conformal factor 
|pOf . One can then write an expansion for to of the form, 

2to = Madm + P 2 L 2 m 1 + ... (8) 

with mi a constant. One can then take the intial data, go{P, L, to) and rewrite it as go{P, L, Madm(P, L, m)) and use 
the above expansion (|8|) for the explicit form of Madm(P, L, to). As a result of this reparameterization, the first and 
second order terms of the initial data are left invariant. That is, one simply takes the initial data and where it read 
"2m" one replaces Madm- For the second order pieces this is also true, the second order pieces of (ra) only contribute 
irrelevant I = terms to second order and do not change the initial data. Summarizing, we construct the perturbative 
initial data, and wherever it said to we replace Madm I % and this is consistent to the order of perturbation theory 
we are considering. Therefore our problem is completely parameterized now by the ADM mass, which also facilitates 
comparison with the full numerical data, which are also parameterized and normalized by the ADM mass. This issue 
was the source of significant confusion in this area initially. In particular, the results of [^) are not properly normalized 
and therefore depart from our predictions in this paper for moderate and large values of the momentum (when Madm 
starts to differ significantly from twice the "bare mass" of the holes). 

The concrete details of computation start with (Q). For R > L and for computations only to second order, one 
needs only the portion linear in L of the extrinsic curvature. Keeping terms only to second order, and taking on the 
right hand side of the Hamiltonian constraint, the form of <j) = 4>bl given in (|J), one gets the 0(L 2 P 2 ) piece of the 
conformal factor, 



by solving, 



6 reg = P 2 L 2 4> {2 > + higher order terms, (9) 



99 , f0 \ / 9 i 9 9 P(l — 2 cos 2 8 2 + 13 cos 4 9) , . 

P 2 L 2 V 2 (2) = S {2) = -72P 2 ! 2 — ; - (10) 

(2P + m) 7 v ; 
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with the boundary condition 4>^ — > at R — > oo. We can simplify the solution of this Poisson equation by 
decomposing the source into multipoles: 



:0) = 


1056 


P 2 L 2 R 
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The solution for the monopole and quadrupole parts are 

,(2)(£=0) _ _ 



11 P 2 L 2 (8Rm + 20R 2 + m 2 ) P 2 L 2 q 



50 (2R + m) 5 R R 

(2)(i=2) _ 8 -P 2 i 2 (™ 4 + 10m 3 i? + 40m 2 R 2 + 80mR 3 + 80i? 4 ) P 2 L 2 q 2 
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(11) 



One can also solve for the I = 4 piece but we will not need it. 

The solution contains two constants qo and q 2 representing the homogeneous solution of the Poisson equation. These 
constants determine, in effect, what boundary conditions are being chosen for the conformal factor. The choice we have 
made is that re gj and hence <^ 2 \ is regular everywhere. The wrong choice of qo or q 2 means that when the solutions 
in ( pi] ) are continued to R < L they will be singular, so that <fi = 4>bl does not contain all the information about 
the singularities. (This would be the case, for instance, if we took qo and q 2 to have the values for the symmetrized 
solution.) 

To determine what values of qo and q 2 give a regular solution, we start by noticing that from (^) we can see that 
the asymptotic form of the conformal factor is, 
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(12) 

On the other hand, we know that the regular part of the conformal factor admits an expansion of the form, 

^ A n (9,P,L) 

- — ep — ( } 

n— 1 

Since the right hand side of ([?]) falls off as 1/ R 6 one can conclude that the first three coefficients of the above expansion 
are part of the homogeneous solution, and have the form 

A (6,P,L)=Q (P,L) 
A 1 (9,P,L) = 

A 2 (8,P,L) = Q 2 {P,L)P 2 {cos{6)) . (14) 

We therefore see that go and q 2 are just the leading coefficients in an expansion of Qo and Q 2 in terms of P and L. 

One can obtain a closed form expression for Qo and Q 2 by applying Gauss' theorem to (^), and using the fact that 
(by choice) <fi Teg is regular on the whole R, 9 plane E 2 . This takes the form 



V0 rcg ■ dS = / S(R, 0, P, L) d 2 S, (15) 
as 2 Jt, 2 

where S(R,9,P,L) represents the right hand side of (|7]), and the integral on the left is evaluated over the boundary 
of the plane at infinity. It is clear that the only term that contributes to this integral is the leading term for the 
expansion of <f>. From there we can therefore determine Qo- Considering the same construction, now for R 2 <j) one can 
determine Q 2 . The results are, 
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(16) 



dr / dO R 4 sm(9)P 2 (cos(9))S(R,9,P,L) 



so therefore for the leading terms we get, 

1 d 4 

QQ = 
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8 8P 2 dL 2 

1 d 4 
8dP 2 dL 2 



d6 R 2 sm(6)P (cos(6))S(R, 9, P, L) 
dd R 4 sm(9)P 2 (cos{9))S{R, 9, P, L) 



(17) 



P=0.L=0 



P=0.L=0 



by. 



These expressions are straightforward to evaluate, especially since to the order of interest we can replace the source 



(18) 



where K 2 is the square of the trace of (||) and <j>o is the conformal factor evaluated for P = 0. The integrals 
however, cannot be solved in closed form. Instead they were computed numerically (in several different ways). The 
numerical treatment of S(R, 9, P, L) requires some care. As pointed out near the end of Sec. IIB, though the source 
term is regular, both the numerator and denominator on the right hand side diverge at the points representing the 
holes. The results we get for the constants are 



qo = 0.219/m 3 



q 2 = 0.224/r 



(19) 



These numbers are in excellent agreement with an approximate calculation due to Brandt and Briigmann pi[ . They 
obtain approximately the correction of the ADM mass due to the momentum in the initial data. The leading term 
in the expansion in L of their formula is precisely our qo. Their result is 11/50 ~ 0.22. One can reproduce these 
formulas by considering expansions of the integrals considered in powers of L. 



D. Casting the initial data in the Zerilli formalism 



Having the initial data for the problem, we now can input it into the perturbation formalism and evolve it. The 
first order perturbations are evolved with a Zerilli equation. The second order perturbations are evolved with a 
Zerilli equation with a "source" term quadratic in first order perturbations. The details of how this is done for the 
momentarily stationary Misner initial data was described in GNPP. Those details, however, were rather specific 
to the Misner case. In particular, the formalism in that work used the fact that in the Misner initial metric the only 
first order perturbations are quadrupolar, and hence the source in the second order Zerilli equation is constructed 
entirely from t = 2 first order perturbations. Those details also assumed that certain of the second order initial metric 
perturbations vanished. 

It will be convenient to use gauge transformations to satisfy these same conditions, so that the previous formalism 
can be used. The initial metric (because it is conformally flat) has the correct Misner-like second order form. The 
extrinsic curvature, however, has a first order t = perturbation which generates £ — perturbations in the evolved 
data. These i = perturbations would contribute to the source term of the second order Zerilli equation. Below we will 
use a first order gauge transformation to eliminate this first order perturbation. This transformation, however, changes 
the second order initial metric, taking it out of "Misner form." We then use a second order gauge transformation to 
restore it to the Misner form. 

Let us start by writing the perturbations in the standard Regge- Wheeler jl5j notation for the multipolar decom- 
position of a metric tensor g a b, ie, 

gu = (1 - 2M/r) ^ H { Q e) (r, t)Pt(cos 9) (20) 

t 

g tr = J2 H i\ r ' t ) p e( co s 6 ) (21) 

i 
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^2h ( i ) (r,t)(d/d0)P i (cos0) (23) 

t 

(1 - 2M/r)- 1 J2H i 2 e \r,t)Pi{cosO) (24) 

i 

r 2 ^ [K W (r,t)+GW(r,t)(d?/d6 2 )\ P e (cos9) (25) 
i 

r 2 sin 2 6j2 \K {l \r,t) +G& (r,t) cot 9(8/86)] P e (cos8) . (26) 



In these expressions, r is related to R, the radial coordinate in the conformally related flat space, by R — (y/r + 
y/r — 2M) 2 /4. The "background mass" M, as previously discussed, is the ADM mass computed numerically for a 
given choice of m, L, P. 

Since the initial geometry is conformally flat, the only non-vanishing perturbations are those in H 2 and K. The 
quadrupole parts, to second order, of these perturbations are: 



H V=*> _ K (i=y _ ^ ^ — " | — (27) 

2 y/f(y/r + \/r-2M) 5 lr(yfr~ + y/r - 2M) 10 yftiyft + Vr - 2Mf 

256[12r 2 - 9rM + M 2 + (8r - 3M) y/r y/r - 2M]L 2 P 2 

35r 3 ( v / r + y/r -2M f 

To describe the perturbations of the extrinsic curvature we shall use a notation like that in (f20|)-(|26|), but shall 
prefix extrinsic curvature quantities with a "if". Thus, for example, K rr — $D<(1 — 2M /^KH^ Pi(cos 6). The 
non- vanishing monopole perturbations of the extrinsic curvature are 

jfflf*) = 2^ (28) 
KK^=-^ (29) 



and the quadrupole perturbations are 



KU ' - ' = — + 16 ( 6r - 11M + G Vr^r-2M) 

2 r3 7r 7/2 +x / r -2M) 5 

KGi ^ ) = _PL + SPL^VF^M 

r3 7r 7/2 (^ +x / r -2M) 



KR(i=2) 5 PL 8 PL 3 (3r-5M + 3^^/7^211) 
7-3 7r 7/2 (^ +v / r _2M) 5 

Ktf* = 32 P£3 (33) 

7y/r-2Mr 3 / 2 (V^+Vr-2M) 4 

We start the process of gauge transformations by writing a general I = and £ = 2 first order gauge transformation 
vector, 

e = M« (fc2) P 2 (cos0) + M« (fc0) (34) 

C = Af«< £ = 2 >P 2 (cos0) + Ml (1)(fc0) (35) 

e , =M ( 1 )(,= 2 )^^) ) (36) 

and we choose all components to vanish except, 

M (D(*=o) = (37) 

MyV - 2M 



M ( 1)(fc0) = _Vr-2M tpL m 
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This gauge transformation eliminates the £ = perturbation of the extrinsic curvature to first order, and leaves the 
I — 2 first order initial data unchanged, but it introduces quadratic changes in the second order components of the 
initial data. To compute these second order changes, we need a four dimensional metric, whereas up to now we have 
only dealt with the initial values. We assume zero perturbative lapse and shift to all orders and use the initial data 
to write an expansion in powers of a fiducial time t for the four dimensional metric around t = 0, 



= {g^)t=o + t{d t g^)t=o + 



(39) 



where (gu V )t=Q is constructed in a straightforward manner with the 3-metric and the chosen lapse and shift, and the 
time derivative of the perturbative piece of the metric h^ v is completely determined by the extrinsic curvature, 



2M 

{dth ab ) t= Q = -2K ab \l 1 . 



(40) 



We then apply the formulas for gauge transformations to the above constructed metric and take the limit t — ► to 
recover the initial data in the new gauge. 

The second order changes due to quadratic combinations of the first order gauge transformation have I = and 
I = 2 components. We will ignore the first, since they are non-radiative. The £ = 2 second order metric that results 
from the gauge transformation of (p^)-(Bq) is: 
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r i M( y /r + Vr-2M) 4 
16PL 3 Wr - 2M 



(41) 
(42) 
(43) 
(44) 
(45) 



In the formalism of GNPP the initial data was taken to have Hq = Hi = ho = up to second order. This was 
true of our perturbed metric before the gauge transformation of (p7|)-(|3q), but is not true of the post-transformation 
metric of (fil[)-(45). We now restore the conditions H = Hi = h a = 0, for the quadrupole, with another, purely 
second order, gauge transformation: 



M, 



(1=2) 



A t 2 PL 2 M ~ 6 L + Pt ^ r ^ 2 M (Vr + v'r- 2Mf 
3 r r Vr - 2M{y/r + y/r^2Mf 



M 



(1=2) 



At P L 2 yfr - 2M 



AL r 3 M - Ptyjr -2M{y/r -2M + y/rf 



r 6 (y/r + ^/r - 2M) 5 



M (^2) = 1 r L 2 M t 3 ~ 8L r3M ^ 1 + Pt(r — 2 M) (y/r + ^~2Mf 
3 r 10 (v^+ y/r-2Mf 



(46) 
(47) 
(48) 



With this transformation, the final form of the first and second order parts of the quadrupole metric perturbations 
read, 



H (i)(.e=a) _ K (i)(t=2) 



16 



ML 



H 



(2)(l=2) 



35 L2P2 



{y/r + y/r-2M f y/¥ 
1272 r 5/2 + 1240 \/r - 2 Mr 2 + 2480 (r - 2 M) r 3/2 + 2480 (r - 2 M) 3/2 r + 



(49) 



+ 1240 (r - 2 M) 2 y/r + 248 (r - 2 M) 5/2 / (y/r + \/r - 2 M^j 5 r c 



M 



+ 



128 L 2 P 2 q2 



192 M 2 L 4 



(y/r + y/r - 2 M)° y/¥ 7 (y/r + \/r - 2 M) 10 r 
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#(2)(*=2) 



642 r 5/2 - 1910 Vr - 2 Mr 2 - 3820 (r — 2 M) r 3/2 - 3820 (r - 2 M) 3/2 r- 



1910 (r - 2 M) 2 V? 7 - 382 (r - 2 M) 5/2 



(v/r + Vr - 2 M^j r 3 M 



128 L 2 P 2 q2 



+ 



192 M 2 L 4 



(V¥+Vr~ 2l) 5 Vr 7 ( V^ + Vr - 2 M) 10 r 
2P 2 L 2 



Mr 3 

and the extrinsic curvature is, 



^(D(^=2) = 



4PL/r 3 
-PL/r 3 

KK (l)(l=2) = _ 5 p L/r 3 



(50) 



Kh\ 



(2)(i=2) _ (108 + 52 - 2 M) PL 3 



iff/. 



(2)(<=2) 



7 (v^ + - 2 M) r 3 / — 2 M 
(296 r + 184 Vr Vr — 2 M + 64 M) PL 3 



KG (2)(£=2) = _ 



7 (V^+ Vr-2M) r 7 /2 
(48 r - 8 V^Vr - 2 M + 16 M) PL 3 



7 (V^+ Vr-2M) r 7 / 2 
(2) _ (4 r + 116 V^Vr — 2 M — 16 M) PL 3 
7 (V^+ Vr-2M) 5 rV2 

For perturbations satisfying the Misncr conditions (Ho = Hi = ho = 0) the first order, quadrupolc, Zcrilli function, 
in the notation, of GNPP is given by 



2M 



3(2r + 3M) 
and its time derivative by, 



rH. 



(l)(i=2) 



+ 3r 2 Gi 1 ^-r 2 Ki 1 ^-6h[ 



(1)(*=2) 



_1_ 1_k-(i)(<=2) 
3 



2 Vr - 2 M 
V^(2r + 3M) 



_ r 2 KK {i)d=2) _ r (3 M _ r) KG (i)(£=2) + (r _ 2 M) A7i< : 



(1)«=2)' 



(51) 



(52) 



Here use has been made of the first order Einstein equations to simplify the occurrence of higher time derivatives. 
With the notation and formalism of GNPP, the second order, t = 2 Zerilli function \ is computed to be 



X 



Ur 4 p {-KG {2) r 2 + r 2 KK {2) + 3rKG (2) M - Kh^ r + 2 Kh { ? ] M 
-2r 9 / 2 (2r + 3M) (X«) 2 - 4r 5 p 3 /4^ (1) + AK^r 6 KK^ p - Ar 4 h ( ^ KK^ Mp + 2r 6 KK^ H { 2 1] p 
+6r 5 p (5r + 9M)KG (1) G (1) + Ar 5 p (-2 r + M) KG {1) H { 2 1] -r 6 p (-3r + 4M) KK^G^ 
+r 6 p (-13 r + 20 M)KG (1) G i J 1) -2r 5 p (8r + 9M) K^KG {1) + 4 

-2r 5 p (7r + 6M)iffi- (1) G (1) +r 5 p 5 ia^G'( 1) + 2rV 3 M? (1) ^ 1} + 4r 5 p 3 ia^ 1) ^ 1) + 4r 5 p 3 ^ 1) ^ 1) 

-6rVG (1) ia« - 2r 7 p 3 KGWKW -2r 7 p 3 KK^G^ -2r 3 p 5 Kh^ r h^ - 2 r 4 p 3 ^ KH { 2 1] 

+r 6 p 3 G^KH { 2 1] +4r 3 p 5 Kh[ 1) h{]l-7r i p 3 Kh{ 1) H ( 2 1) +2p 3 r 2 (7 M + 22 r) Kh^h^ 

-8 x A : (2r + 3M)/) 4 (/ i < 1) ) 2 +rV (M — 26 r) Kh^G^ - 2A r 5 p 3 KG^ h^ + 12 r 5 p 3 h { ^ r KG^ 

+ 12 r 7 p 3 KG^G^ + 2 r 6 p 3 KG^H^ - 4 r 9 / 2 (2 r + 3 M) p 2 G«X« + 8 r 5 ' 2 (2 r + 3 M) p^G^ 

+8 r 5 / 2 (2 r + 3 M) p 2 h { ^K^ - 2 r 9 / 2 (2 r + 3 M) p 4 (G^) 2 ] [7 r 9 ' 2 (2 r + 3 M)]" 1 



(53) 
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where p — yr — 2Af, and the time derivative of the second order, I = 2, Zerilli function is given by, 

Xit = -4r 15/2 {7 M -6r) p^H^G^ + 16 Mr 7 p 5 (2r + 3 Af) K {1) ffG W - 4r 17/2 p 10 KG^ KhQ (54) 
+28 H ( 2 2) r 13/2 (Af + r) p 6 + 2 r 13/2 (18 Af 2 + 17 rM - 8 r 2 ) p 6 G^K^ + 4 r 13/2 (5 r + 9 AT) p 6 G [1) H { 2 1] 
-4 r 15 / 2 (M + 12 r) p s KG (1) Kh ( ^l + 4 r 13 / 2 (2 r + 3 Af ) p s G {1) H ( 2 1] r + 4 r 15 / 2 (3 M + 50 r) p 8 KG^Kh^ 
_4 r i7/2 (17 m+ 15 r) p 6 (fiG (1) ) 2 + 8r 19/2 (3M-r) p 6 KG w KK ( r 1) + 4r 13/2 (6M-5r) p 6 G (1) ff (1) 
-12 r 11 / 2 (2 r + 3 M) p^^W _ 48 r 21 / 2 p 8 (KG^) 2 - 116 r 13 / 2 p 8 (ift^ ) 2 

-28 r 15 / 2 p 6 f^ 2 > - 12 r 19 / 2 ^ 6 ^^)) 2 + 14 r^' 2 p 6 {K^f - 6r 13 ' 2 (12 Af 2 + 19rAf - 6r 2 ) p e G^G^ 

-24 r 9 / 2 (10M-3r) p 8 /^/^ + 4 r 19 / 2 (2M-3r) p 6 KG^KK^ + Ur 15 / 2 Mp 6 K^ 

-12 r 7/2 (42 M 2 - 14 rM - 3 r 2 ) p e (h[ 1] ) 2 + 3 r 15/2 (54 M 2 - 30 rAf + r 2 ) p e (G { r 1) ) 2 

+2r 17 ' 2 p*H 2 [ ] ) r K^ - 28 r 9 / 2 (6M 2 + 5rAf - 2r 2 ) p 6 h 2 - 3 r 17 / 2 p 8 (ff «) 2 

_ r i3/2 (2 M — 5 r) p 6 (fff >) 2 - 28 r 11 ' 2 (2 r + 3 M) p 8 h^ 2 + 4 r 11 / 2 (4 Af - 7 r) p'h^H^ 

+16 r 21 / 2 p 8 KG^KKp +48 r 15 / 2 p 6 (G«) 2 + 16 »V (2r + 3 Af ) KK^K^ 

+64 r 4 p 9 (2 r + 3 Af) K7i^ /i^ + 16 r 8 p 7 (2 r + 3 Af ) KK {1) G^ - 32 Af r 5 p 7 (2 r + 3 M) h[ 1] 

+16 r 8 p 7 (2 r + 3 A/) KG^K^ + 2 r 15 / 2 (3 Af - r) p^G^H^l + 16 r 7 Afp 7 (2 r + 3 Af) KG {1) G { ^ 

-32 r 6 p 9 (2r + 3 M) Kh^G^ - 32 r 6 p 7 (2r + 3Af) KK^h^ + 16 r 8 p 9 (2 r + 3 M ) KG^ G W 

-32 r 6 p 7 (2 r + 3 Af) Kh^K^ - 12 r 13 / 2 (2 r + 3 A/) (Af - r) p 6 G^K^ - 40 r 11 ' 2 p s h^ ] 

+12 r 13 / 2 (10 Af - 3 r) p 8 G« Tig + 12 r 11 ' 2 (8 Af 2 + 11 rAf - 7 r 2 ) AiM* 5 

+8 Af r g / 2 (AM - 7 r) p 6 h^H^ + 36 r 9 / 2 (4 A/ 2 - 2 r 2 + rAf) p^G™ 

+4r 13 / 2 (5 A/ 2 + 2 rAf - 2r 2 ) p e KG^Kh^ ] - Ar 11 ' 2 (3 M - r) p* H^h^ - Wr 19 / 2 /) 8 ^^ 
-4r n / 2 (35 M 2 - 19 rAf + r 2 ) p 6 K^h^ ] + 16 r 17 / 2 p s Kh^ r KKW _ 2r 15 / 2 (27Af 2 - 7r 2 ) p 6 G«f^ 
+8r 13 / 2 p 8 / l Wff + 32 r u ^ 2 p s Kh^KH^ - Ar 13 ' 2 (10 Af - 3r) Ai>#« 
+16 r 15/2 Mp 6 Kh^ ) KK w +4r 17/2 (Af + 12 r) p 6 KG w KH 2 [1) 

-4 r 15/2 p 6 ^ (1) ff^ 1} - 14 r 13/2 (6 Af 2 - 3 rAf - 2 r 2 ) p 6 G< 2) + 4 r 17/2 (17 Af + 16 r) p 6 KG (1) KK w 

-~Ar 15 ^ 2 Mp e K^K^ -8r 13 / 2 p 10 Kh ( llKh ( ^ - 2r 15 / a p 8 iT 2 (1) fl£> 

_ 4r i9/2 (_i 3r + 20 M )p e KG^KG^ -%r w l 2 p 6 KH^KK^ - 32 r 17 ' ' 2 p s KK^ Kh{ 1} 

-S2r 6 p 9 {2r + 3M)KGi 1) h ( i L) ] [l4 p 4 (2 r + 3 Af ) r 17/2 j 1 

where a subscript r denotes differentiation. To arrive at the expressions in ([53] ) and (|54|) the second order Einstein 
equations have been used to eliminate higher order time derivatives. The above expressions were automatically 
computed with Maple computer algebra codes. It is impractical to give more details of their construction in print. 
The source codes and documentation can be found in our anonymous ftp server Jig ], 

When the explicit 3-geometry and extrinsic curvature of (|49|)-(po|) are put into the expressions of (jfij) and ([54]), we 
arrive at the following initial data for the first and second order Zerilli equations, 

8 ML 2 (5Vr - 2Af + 7^/r) r t 

ipt=o=»- , 5 (55) 

3 (Vr + Vr - 2Af ) (2r + 3Af ) 

s/r - 2MPL(8r + 6M) 

= rS/ 2 (2r + 3Af) (56) 

_ 512 Af 2 f 4 + 16 (9y^ + 17vV - 2M ) Vr - 2MMPL 3 

7 (v^+Vr-2M) 10 r 7 {JT+^r- 2M f (2r + 3Af )r 5 / 2 
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64 M 2 (10 r - 10 M + 38 ^r~y/r - 2 M) y/r-2 ML 4 

~Kt— o = — 

7 (Vr+ V?--2Af) 10 (4r + 6M)r 5 / 2 

64 (4r + 14M)Vr- 2M M P L 3 Q ^M (py/r — 3\/r — 2M) \Jr — 2M L 2 P 2 q2 
7 (y/¥+Vr-2M) 5 r 5 (y/f + Vr - 2 uf (2r + 3 M ) r 5 /2 

16 

+— Vr-2ML 2 P 2 (1750M 4 - 9849rM 3 + 2331r 2 Af 2 + 7182r 3 Af - 2892r 4 
35 v 

-Vta/t - 2Af(3148r 3 - 4130r 2 Af - 4935rAf 2 + 4375A/ 3 )) /((2r + 3M) (y/r + Vr - 2 A?) 5 r 6 ^ . (58) 
We are now ready to evolve the initial data and compute waveforms and radiated powers. 



III. EVOLUTION 



A. The Zerilli equations 



The initial data generated in the previous section is now fed to the first and second order Zerilli equations 

where r* is the usual "tortoise" coordinate covering the exterior of the black hole, 

r* = r + 2A/ln(r/2A/- 1) (61) 

so the horizon is at = — oo and spatial infinity at r, = oo, and the potential and source terms in the Zerilli || 
equations are given by, 



V(r) 
S 



= 1 



2AA f4r 



- 2 r 72M 3 12M„ w , J 3M 
-5-(2-l)(i + 2) 1- — 



12 ^ 3 



A 2 



12(r 2 + Mr + M 2 ) 2 



2(Z-l)(i + 2)1(1 + 1) 



W»,t) -4 



(2r 3 + 4r 2 A/ + 9rAf 2 + 6M 3 ) 



rV 2 A r 6 A 
(112r 5 + 480r 4 M + 692r 3 M 2 + 762r 2 M 3 + 441rAf 4 + 144M 5 ) 

r 5 /i 2 A 3 

18r 3 - 4r 2 M - 33rM 2 - 48A/ 3 



rA 

1 



(62) 



3r 4 ^ 2 A 



, , 12r 3 + 36r^M + 59rAP + 90A/ 3 . , , 2 

VV + 3^ ) 



(2r 5 + 9r 4 M + 6r 3 M 2 -2r 2 Af 3 -15rM 4 - 15M 5 ) , 2 (r 2 + rM + M 2 ) , , 
hl2^ 5-s-: — V,t-0,tr 



r 8 fi 2 A 



r 3 fi 2 



(32r 5 + 88r 4 M + 296r 3 Af 2 + 510r 2 Af 3 + 561rAf 4 + 270M 5 ) , , 1 , , 

~ 2 r7 ^ A2 Wir + 3^2 Vir ^,trr 



2r 2 - Af 2 
r 3 yuA 



3r 2 + 12rM + 7AP 



r 4 /iA 



4(3r 2 + 5rA-f + 6M 2 ) 
37 s 



^iA 



3r - 7M 
3r 3 /Lt 
2r + 3M 



-ip, r ip,t 



M 
r 3 A 



(63) 



where A = (Ar + 3M) with A = (£-!)(£ + 2)/2 and fi= (r- 2M). The potential is given for general ^ in @ but 
we will use it only for I = 2. 

As can be seen, equations (^) and (^) have the same form, including the same potentials, but the second order 
equation has a source term that is quadratic in the first order Zerilli || function and its time derivatives. We have 
written a fortran code to evolve these equations by a simple leapfrog algorithm. Convergence to second order was 
checked and special care was taken to avoid noise from the high derivative order of the source term. 
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To find the gravitational waveforms and power a transformation must be made to a gauge that is asymptotically 
flat to first and second order. The details of this process were discussed in GNPP and will not be repeated here. 
The result is that the transverse-traceless perturbations, in the asymptotically flat gauge, correct to second order, are 
encoded in the quantity 



ft 4-\ 9 ^ 



Id f dip 



(64) 



(where it is understood that all quantities are I = 2) and this is the quantity we shall plot below when we give 
waveforms of the outgoing gravitational radiation. The first order part of the radiation is given by the leading term in 
(|64|); the terms in square brackets are second order. From the Landau-Lifschitz pseudo-tensor in the asymptotically 
flat gauge (as discussed in GNPP) we find that the radiated power is 



p 3 ffy 

Power = To 1 at 



Id ( dip 



(65) 



(Note that the perturbation parameter e that appeared in |Q is now incorporated into the definition of the Zerilli 
functions we have used in the paper, as can be seen in formulas (|55|-j5"8|) . We have also directly computed the 
"renormalized" second order Zerilli function in ([53])). 

Before we move on to present our results and compare with the numerical relativity simulations of the Pots- 
dam/NCSA/WashU group (see BAABRPS), it is worth pointing out, again, that the numerical relativity simulations 
are for "symmetrized" initial data, in which an infinite number of "image charges" is used to construct initial data 
representing two throats connecting two isometric asymptotically flat universes. By contrast, the problem we are 
solving corresponds to three asymptotically flat universes. In the limit of zero momentum the numerical simulations 
correspond to the 1960 Misner J?J initial data, and our results correspond to the Brill-Lindquist JHJ initial data. 

For the range of separations we are going to discuss the discrepancies between these two types of data are insignifi- 
cant. (Although we are working in the "close-limit," we will consider sets of data far apart enough to make the extra 
terms arising from symmetrization very small), but since the problem is a multi-parametric one, it is not obvious that 
this is true in all the ranges of parameters we will be discussing. More careful studies will be needed if one wants 
higher accuracies than the ones we are going to discuss here. We have also modified the Potsdam/NCSA/WashU 
code to run for unsymmetrized data, and for limited tests the results agree very well with the symmetrized ones in 
the range we are considering. This situation arose due to historical reasons: the numerical code was written before 
our work with non-symmetrized boundary conditions, whereas perturbation theory becomes very cumbersome if one 
starts carrying around the extra terms due to symmetrization. 

One particular problem that one faces when comparing Brill-Lindquist (unsymmetrized) and symmetrized data 
sets is that the sets are parameterized in different ways. There is therefore ambiguity in how to compare the results. 
Abrahams and Price [ fl7| have discussed this in some detail, and show that there are different identifications one can 
take that yield sensible results along a good range of parameters, so we will not repeat the discussion here. We just 
state the convention we are following: For one of our results with momentum parameter P and throat separation 
L, we compare a numerical relativity result with the same ADM mass and same numerical value of the momentum 
parameter, and with a separation parameter given by 

L = 2v/4K 2 (/io) ■ (66) 

Here /xo is a parameter originally introduced by Misner |Q that is commonly used to parameterize symmetrized binary 
black hole initial data sets, and 

1 ^ (cothn^o)^ tR „, 
(4Si) + ^ sinhn/io 

oo 1 

s i^E-t — • ( 68 ) 

smnn/zo 

With these choices, the radiated waveforms agree very well when P = 0. Notice that the discussion of Abrahams and 
Price JT7| is only for the P = case. The "best" identification between symmetrized and unsymmetrized data could 
probably be a P-dependent notion. We will ignore this issue here, but it clearly requires further study. 
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B. Fixing t = 



In the formalism of GNPP we chose to fix the coordinates by requiring that the metric be in the Regge- Wheeler 
jTJ| gauge to first and second order. This can always be done, but it turns out that the coordinates are not quite 
uniquely fixed. The problem is quite generic and it has to do with how perturbation theory handles time transla- 
tions in situations where the background spacetime is time-translation invariant. Consider an exact quantity f(r,t) 
approximated by a perturbative series expansion, 

f(r, t) = /(°) (r) + e/« (r, t) + e 2 / (2) (r, *) + ... , (69) 

and perform now a first order gauge transformation corresponding to a pure time translation, 

t -> t' = t + ec, (70) 

with c a constant, independent of t and r. Replacing if by the above expression (and noticing that dt = dt'), we get 

/M') = / (0) M + e/ (1) M') + 6 2 / (2) M') + 0(e 3 ) 

= / (0) (r) +e/ (1) (r,t) + e 2 (c/«(r,i) + / (2) (r,<)) + 0(e 3 ). (71) 

So we see that the "second order term" in the expansion of the metric depends on the origin chosen for time. If one 
starts with perturbations in the Regge- Wheeler gauge, a transformation of type ( |70"| ) leaves the perturbations in the 
Regge- Wheeler gauge, but the second-order metric is changed, and in fact depends on an arbitrary constant c. This 
indicates that a comparison of quantities to second order in perturbation theory around stationary backgrounds can 
be quite misleading: the same metric can have very different second order terms depending on the origin of time 
chosen. Worse, these terms can be quite large, and are completely artificial. 

It is interesting to notice that if one computes the radiated energies using the formula we discussed previously (|6^ ) . 
the results are unchanged — as expected — by time translations (the additional term turns out to be a total derivative 
that does not affect the computation of energies). But we want to go beyond giving perturbative results for radiated 
energy. We want also to compare perturbative waveforms with those of numerical relativity. Since these waveforms 
are second-order correct quantities given as a function of time at a particular "observation" radius, we must be sure 
that we are using the same zero of time for both waveforms, that from perturbation theory and that computed with 
numerical relativity. 

Fortunately, it is not difficult to eliminate the time-shift ambiguity in the metric. To do this we separate the 
waveform f(r,t) given in ( |6~4| ) into first and second order parts f^- 1 ' and and we construct the quantity 

r^/WM);( 2 )M) 

co = 72 — ■ ('2) 

p oo dt[fW(r,t)_ 

We then perform the time translation t — > t + cq, arriving at the "physical" value for the second order waveform 

f$ ys (r,t) = fW{r,t)-c fW(r,t) . (73) 
Equivalently, we adjust the zero of time, and hence f( 2 '(r,t) until the integral in the numerator of (173) vanishes. The 

— (2) 

same coordinate fixing must be done to the numerically computed waveform f num . To do this we define /n Um to be 
/num — f^- We then adjust the zero of time so that the integral of /num/^ 1 -* vanishes. 

These observations about time-shifts are also true in the time-symmetric case. We have recomputed the results of 
Pj with the zero of time fixed as above and have found that the results are changed by less than 1%. For boosted 
black holes, on the other hand, this time fixing is crucial for seeing the high accuracy agreement of the perturbative 
and numerical relativity results. 



C. Results: radiated energies 

We start to summarize our results by computing the radiated energy as a function of momentum for head-on 
collisions of black holes released from a separation of fio = 1-5, L/(0.5Madm) — 4.2. The results are depicted in 
Fig.| 
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The figure shows the characteristic "dip" at low values of the momentum that was first noticed in BAABRPS. 
An important difference between that paper and the present results is that here, as explained in Sec. 11C, we are 
normalizing both the numerical and the perturbative results using the same ADM mass. This leads to a much better 
agreement for large values of the momentum than that observed in BAABRPS. As an example of the size of the 
difference, for P/M ADM = 1, P/(2m) « 3 and for P/M ADM = 3, P/(2m) w 15. 

A remarkable fact is that first order perturbation theory agrees very well even for large values of the momentum, and 
second order perturbation theory confirms this fact. This at first seems puzzling since our initial data was obtained 
through a "slow" approximation in which the momentum was assumed to be small. However, as was observed in 
BAABRPS, for large values of the momentum the initial data is "momentum dominated", meaning that the extrinsic 
curvature completely dominates the initial data. Therefore the errors made in computing the conformal factor via 
the slow approximation become less relevant than might be supposed. 




FIG. 1. Radiated energy in head-on black hole collisions as function of the momentum for a separation of fio — 1.5, 
L/(0.5Madm) = 4.2. Depicted are the close-slow approximation and the full numerical results of the Potsdam/NCSA/WashU 
group. Even for large values of the momentum, the first order results overshoot and the first plus second order undershoot the 
numerical results by only 20%. The inset shows the "dip" region. 

The overall picture of the energy therefore is very encouraging, the approximations presented seem to be working 
even beyond their expected realm, and second order perturbation theory is capable of tracking this fact, playing the 
expected role of "error bars." This approach is not without pitfalls, however. In order to illustrate these, we turn 
to Fig.||, which shows a close up look at the energy picture and also includes results for black holes initially boosted 
away from each other. 

The first thing we notice is, that for black holes boosted away from each other there is, as expected, no "dip" in the 
energy. The dip is a first order effect that is due to a cancellation between terms that are momentum independent and 
terms that are linear in momentum. The cancellation turns to addition in the case of negative (outwards) P. We also 
see that first order calculations are less accurate at the dip than at higher values of the momentum. This is somewhat 
puzzling since our approximation should work better the smaller the momentum. What seems to be happening is 
that first order theory does not accurately reproduce the higher order terms that make important contributions to 
the energy after the leading terms cancel in order to produce the dip. This is confirmed by the fact that first plus 
second order results are indeed very accurate at the dip. 

An instructive feature of these results is that for black holes boosted away from each other a cancellation of the 
second order terms takes place around P/(2m) = 0.9, P/Madm = 0.36. Clearly one cannot regard second order 
perturbation theory as giving error bars when it is cancelling out. Moreover, it shows that second order results beyond 
that value of P can only be taken as rough indicators. We will return to this cancellation in somewhat more detail in 
connection with waveforms. 

Another issue to be mentioned is how crucial it is to have chosen the ADM mass of the initial slice as the mass 
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0.0 1.0 

PM WM 

FIG. 2. We depict energies for black hole collisions in which black holes were initially boosted towards each other and also 
for the case in which they were boosted away from each other. The "dip" effect is only present when the black holes are moving 
towards each other. The approximation given by first order perturbation theory is slightly worse in that case, since leading 
terms in the calculation are cancelling each other to produce the suppression. We also see the cancellation of second order 
contributions that takes place for the outward pointing momentum case. There, first order results undershoot the results for 
small values of P and overshoot it later, the second order corrections vanishing at the point of crossing. 

of the background spacetime used in the perturbative calculations. Our previous (first-order) work on boosted black 
holes used the "bare" mass (ADM mass for P = 0) for the background. This is quite visible if one compares Fig. |l| 
with Fig. 2 of BAABRPS. In the latter, first order perturbation results appeared to disagree with the numerical results 
by over an order of magnitude for P/Madm = 3. That was entirely due to the poor choice of background mass. In 
the present paper, using the numerically computed ADM mass of the initial data, we see that first (and second) order 
results differ by only 20% from the numerical results at P/Madm = 3. 

D. Results: waveforms 

Let us now turn to the examination of waveforms. The numerical code of the Potsdam/NCSA/WashU group 
extracts waveforms at slightly different values of the radial variable for varying P's. We took this effect into account 
and extracted perturbative waveforms at the same radii as was used for the numerical relativity work. In all cases the 
full numerical code has a very limited range of spacetime covered in the evolution. This forces the extraction to be 
done in a rather small range of radii around 20Madm or so. With perturbation calculations we could have extracted 
much further away, but we performed the extraction at exactly the same radius as those used by the numerical code. 
Waveforms were observed to change shape rather significantly from one extraction radius to another even in such 
a close range, but we observed that as long as we extracted the perturbative waveform at the same radius as the 
full numerical result (as opposed to, say, extracting farther out and then shifting the result back) the agreement was 
roughly independent of extraction radius. However, this starts to hint at a main problem in comparing waveforms: 
one needs not only to match amplitudes but it is also crucial to match the phase, at least if one is interested in high 
accuracy. The phase is determined by, among other things, the extraction radius. Determining the extraction radius, 
in turn, requires knowing the ADM mass (since one measures radii in units of ADM mass). Our full numerical code 
for computing the ADM mass, in its present implementation, is accurate to a few percent. (This could be made better 
with more computer power than what is presently available to us; the runs we made had 300 radial zones and 30 
angular zones.) This limits the accuracy with which we know the ADM mass, and hence the accuracy with which 
we can determine the phases. The technique, discussed in Sec. IIIB, of fixing the zero of time is helpful in giving an 
objective way of comparing phases. 
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Let us turn to the results. We present, below, the results for the waveform /(r, t) as defined in (|6J). This is directly 
comparable (up to a time derivative) with the output of the full numerical relativity code, which outputs a Zerilli |5j 
function via the radiation extraction technique of assuming that the spacetime is a perturbation of Schwarzschild and 
reading off the perturbations from the full numerical results. 




FIG. 3. Comparison of waveforms for large values of the momentum, P/(2m) — 14, P/Madm = 2.44. There is good 
overall agreement, but there is some slight disagreement in the details of the waveforms. As one expects for large values of 
the parameters, second order perturbations can at most be regarded as an estimator of error, rather than a way to improve 
the accuracy of the waveforms. It is still remarkable that perturbation theory would work well for such a large value of the 
parameter P. 

Our presentation of waveform comparisons starts with the most disfavorable cases and moves to more favorable 
ones. Figure || shows the comparison of waveforms for P/(2m) — 15, P/Madm ~ 2.44 and Fig.[| corresponds to 
P/(2m) — 5, P/Madm ~ 1.32. As we see, there is very good overall agreement. Notice that (taking into account the 
"time-shift" gauge fixing discussed above) our procedure in the end has no free parameter, i.e., phases and amplitudes 
are predetermined in all cases, which makes the agreement more remarkable. If one looks carefully at the curves in the 
inset, which enlarges the region around the second positive peak, one sees that there are slight phase and amplitude 
disagreements. First order results tends to overshoot the waveform, whereas adding the second order correction tends 
to undershoot. There are slight differences in shapes as well. For large values of the momentum, we can take second 
order predictions as "error bars" only. However, for intermediate values, it is quite clear that first plus second order 
calculations offer a very accurate prediction of the waveforms. 

The reader should exercise care when comparing the results for waveforms with those of energies. This is due to 
a peculiarity of the formula for the radiated power (|65j). As discussed in GNPP, the square that appears in (|6^) 
involves terms that are of "third order" in perturbation theory. Therefore, to keep things consistent, when squaring 
the expression in curly braces, we only keep the mixed term and omit the term that is the square of the second order 
part of /. As a consequence, the second order correction for the radiated energy depends mostly on correlations of 
phases of the first and second order waveforms rather than on their amplitudes. For instance, for the case we are 
studying P/(2m) — 5, P/Madm ~ 1.32, the second order waveforms are only slightly smaller than the numerical 
ones, but the computed energy is 12% lower. 

We now turn our attention to the area of the dip, P/Madm ~ 0.05, P/(2m) = 0.12. In Fig. [| we show the 
waveforms for the inward boosted case (the case with a dip in the energy). We see that second order corrections 
improve the accuracy markedly. Clearly there are strange effects taking place for this value of the parameter. In 
particular, it should be noticed how first order theory overshoots the waveforms rather significantly in the second 
and third positive peak of the waveform, but not in the first one. In view of the fact that the energy is given by the 
correlation of the first and second order waveforms, those discrepancies in the first order waveforms would seem to 
be responsible for the large relative error in the calculation for the energy, even to second order. This is so, in spite 
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FIG. 4. Comparison of waveforms for intermediate values of the momentum, P/Madm ~ 1.32. Here one can see that second 
order theory not only acts as an "error bar," but when added to the first order calculation, actually allows a reasonably accurate 
prediction of the waveforms. 

of the fact that second order calculations yield very accurate waveforms. Figure || shows the case of Pj (2m) = —0.9 
(holes moving initially apart). As could be predicted from the energy plot, a cancellation of the second order terms 
is taking place. In this case, therefore, one cannot regard second order corrections as "error bars," since it is clear 
that higher order terms are important. It is worthwhile pointing out that the cancellation is highly nontrivial, the 
initial data having the same amplitude for both inward and outward momenta. The cancellation takes place in the 
evolution, with the source terms of the second order Zerilli [pj equation playing a significant role. A simple way to 
understand the cancellation is to break up the evolution into three separate Zerilli lq] equations with three different 
initial sources, proportional to L 4 , PL 3 , and P 2 L 2 respectively. What one sees is that the cancellation occurs between 
the PL 3 term and the other two, and clearly depends on the sign of P (for our simulations negative P is outward 
pointing). One can then infer that there is a curve of cancellations in the P,L parameter space that isolates a region 
in parameter space where second order perturbation theory does not help. One cannot reach points in that region 
unless one changes the relative counting of powers of L and P in perturbation theory. A further study of this issue 
could therefore yield interesting results. 



We have seen that the use of combined first and second order perturbation theory can give excellent results for 
waveforms and energies of radiation emitted in the head on collision of two equal mass, initially boosted, black 
holes. The results show, however, that there are some subtleties, not previously appreciated, in the use of higher 
order perturbation theory and in the comparison with results from numerical relativity. The following points deserve 
attention, especially in connection with the application of higher order perturbation theory to further problems. 

a) The comparison of perturbation results and numerical relativity results has pitfalls when comparisons are made 
between problems that are not identical. In our case we compared our perturbation result for an "unsymmetrized" 
(Brill-Lindquist |[o| type) initial data, with numerical relativity results for "symmetrized" (Misner Q type) initial 
data. Had we been comparing with unsymmetrized initial data, the parameters m, L, P for the data sets would have 
had identical meaning. Since the data sets were not identical, a mapping of one parameter set to the other had to be 
imposed. One degree of freedom in this mapping was subsumed in the choice to compare cases of equal ADM mass, 
but the remaining element of choice in the mapping is a source of uncertainty in the high accuracy comparisons we 
are making. (We emphasize that the choice of mapping was made before any results were considered; there was no 
"fine tuning" to improve the comparison. The excellent agreement between the numerical and perturbative results 
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FIG. 5. Waveforms in the "dip" region, P/Madm ~ 0.05, for inward boosted black holes (the case where there is a "dip" in 
the energy). Second order corrections improve the accuracy in the amplitudes significantly, but first order calculations exhibit 
an erratic behavior. Since the energy arises as a correlation between both waveforms, this translates itself into a rather large 
relative error in the perturbative computation of the energy. Notice that waveforms for this particular value of the momentum 
differ markedly from waveforms for other values, the "enveloping" amplitude curve decreasing much more slowly, allowing 
several oscillations to be visible. 




FIG. 6. Waveforms for P/(2m) = —0.9, P/Madm ~ —0.36, i.e., for outward boosted black holes (the case where there is no 
"dip" in the energy). Second order corrections cancel out and therefore are not reliable as "error estimators" nor to improve 
the accuracy of first order calculations. (The first and first plus second order curves are both plotted but are indistinguishable). 
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then must be considered to be, among other things, an indication that there is no great sensitivity to the manner in 
which this mapping of parameters is done.) 

b) There is no unique result that is correct to second order. Different ways in which details are handled will produce 
results that are the same to second order, but differ at higher order. These different results can have different ranges 
of validity and can exhibit different accuracy when compared with numerical work near the limit of validity. One 
example of this feature of perturbation theory is the dependence on parameterization [fl7f . In our perturbative results 
we have seen another simple example: The second-order correct waveform consists of a first order and a second order 
piece. When radiated energy is computed by squaring this waveform one can choose simply to take the square, or to 
truncate the result and omit the fourth order contribution arising from the square of the second order contribution 
to the waveform. (We have made the latter "conservative" choice.) Both results, of course, are equally justifiable for 
the order of perturbation theory we are doing, but the results are noticeably different. 

c) In the present paper we have seen a particularly interesting example of the importance of higher order terms 
and the detailed way in which perturbation theory is applied. To make the comparison between symmetrized and 
unsymmetrized initial data we found that it is important to compare cases of equal ADM mass, but the ADM mass (for 
fixed m) varies quickly with increasing initial momentum. If one computes this momentum dependence perturbatively 
the agreement of perturbation theory and numerical relativity is limited. With the ADM mass computed exactly 
(i.e., numerically) the agreement is greatly improved. This suggests that an a priori physical understanding of the 
dependence on the perturbations can be a very useful guide to an efficient perturbation scheme. 

d) In addition to the numerical computation of ADM mass, another useful new technical detail was developed in 
the present work. A method was found of fixing the zero of time in the same manner for both perturbative and 
numerical waveforms. This fixing of zero had not been important in previous perturbation studies, but was crucial to 
comparison of waveforms for initially boosted holes. 

e) Perturbation analysis in the present paper was carried out for both small separation and small momentum ( "the 
close slow limit"). This makes it particularly difficult to unravel the sources of disagreement with numerical results 
when anomalous cancellations (like the "dip") occur. Although perturbation results end up in excellent agreement 
with numerical relativity results, a perturbative analysis based on small L, but without small P (especially if it could 
be compared with numerical results for unsymmetrized data) , might be useful in improving our understanding of the 
nature of errors. 

f) The current state of the art of numerical relativity presents limitations, both in accuracy and in range of 
simulations of the codes. As a consequence, we were limited to comparing waveforms which are not really in the 
radiation zone. This is a dangerous exercise when it comes to second order perturbation theory. In particular, the 
formula for the radiated power (from which we extracted the concept of second order waveform) assumes that one is 
in the radiation zone. This is true also of the extraction techniques used in the numerical codes to produce a Zcrilli 
H function as output. In short: with the current limitations we cannot rule out that the discrepancies we see in 
waveforms and energies might be within the error margins of the numerical results. 

A general conclusion of this work is that the synergy between numerical results and perturbative calculations will 
probably be one of the major tools that we will have to use to address with any accuracy the problem of the collision 
of two black holes in general relativity. We see this taking place right now. 
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